home *** CD-ROM | disk | FTP | other *** search
/ MacHack 2000 / MacHack 2000.toast / pc / The Hacks / MacHacksBug / Python 1.5.2c1 / Extensions / Numerical / Demo / sieve.py < prev   
Encoding:
Python Source  |  2000-06-23  |  2.0 KB  |  77 lines
  1. #! /usr/bin/env python
  2. from Numeric import *
  3.  
  4. def sieve(max):
  5.         numbers=range(max+1)
  6.         size=int(math.sqrt(max))
  7.         if size<5:
  8.                 trials=[2,3]
  9.         else:
  10.                 trials=sieve(size)
  11.         for i in trials:
  12.                 try:
  13.                         j=i*i
  14.                         while 1:
  15.                                 numbers[j]=0
  16.                                 j=j+i
  17.                 except IndexError:
  18.                         pass
  19.         return filter(lambda x:x>1,numbers)
  20.  
  21. def factor(num,upto=1000000):
  22.         for p in sieve(upto):
  23.                 if num%p==0:
  24.                         print "Can divide by %d"%p
  25.  
  26. def asieve(max,atype='l'):
  27.         times = []    
  28.     #times.append(time.time())
  29.         size=int(sqrt(max))
  30.         if size<5:
  31.                 known_primes=[2,3]
  32.         else:
  33.                 known_primes=asieve(size)
  34.  
  35.     #Initially assume all numbers are prime
  36.     numbers=ones([max+1], 'b') #atype)
  37.     #times.append(time.time())
  38.     #add(numbers, array(1, 'b'), numbers)
  39.     #times.append(time.time())
  40.     #0 and 1 are not prime
  41.     numbers[0:2]=0
  42.     #print trials
  43.         for i in known_primes:
  44.             #Set multiples of i to be nonprime
  45.             numbers[(i*i)::i] = 0
  46.         #times.append(time.time())
  47.     #Those entries which are nonzero are prime
  48.     #a = add.accumulate(ones([len(numbers)]))-1
  49.     a = arange(len(numbers))
  50.         #times.append(time.time())
  51.         r = repeat(a, numbers)#numbers)nonzero(numbers)
  52.         #times.append(time.time())
  53.     #times = array(times)
  54.     #print times[1:]-times[:-1]
  55.     return r
  56.  
  57. def afactor(num,upto=1000000):
  58.         #Calculate all of the primes up to upto
  59.         s = asieve(upto)
  60.         #Return just those primes that divide evenly into num
  61.         for d in take(s, nonzero(equal(num % s, 0))):
  62.             print "Can divide by %d"%d
  63.  
  64.  
  65. import time
  66.  
  67. #start = time.time()
  68. #factor(129753308,upto=99999)
  69. #stop = time.time()
  70. #print "factor took %7.3f seconds"%((stop-start))
  71. start = time.time()
  72. #asieve(999999)
  73. afactor(129753308,upto=999999)
  74. stop = time.time()
  75. print "afactor took %7.3f seconds"%((stop-start))
  76.  
  77.